/*******************************************************************************

This code file defines bootstrap utilities for the marginal fiscal cost 
and elasticity

*******************************************************************************/

xtset,clear

*** Average fiscal cost

	capture program drop boot_acost
	program define boot_acost, rclass
				 		  
		preserve
								
		logit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict xb, xb
		
		gen pr = exp(xb)/(1+exp(xb))
				
		gen mvalue_bldg = (assesstot - assessland)/(frac_assess*underassess)
			
		gen assessprop = pr * mvalue_bldg * dtaxrate_onsite
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		
		collapse (sum) assessprop inclusionaryunits
		
		gen acostperunit = assessprop/inclusionaryunits
		
		summ acostperunit
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Average fiscal cost, reweighted

	capture program drop boot_acost_rw
	program define boot_acost_rw, rclass
				 		  
		preserve
								
		logit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict xb, xb
		
		gen pr = exp(xb)/(1+exp(xb))
		
		gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
		
		gen mvalue_bldg = (assesstot - assessland)/(frac_assess*underassess)
			
		gen assessprop = pr * mvalue_bldg * dtaxrate_onsite
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		
		collapse (sum) assessprop inclusionaryunits [aw=reweight]
		
		gen acostperunit = assessprop/inclusionaryunits
		
		summ acostperunit
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Average fiscal cost, GMM

	capture program drop boot_acost_gmm
	program define boot_acost_gmm, rclass
				 		  
		preserve
				
		logit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b_init = e(b)
				
		gmm (inclusionary_onsite - invlogit({xb:dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit constant})), ///
			instruments(dtaxrate_onsite_sim $controls_lot_bs $controls_block i.borough i.yearpermit) vce(cluster nta_boot) ///
			conv_maxiter(10) from(b_init) quickderivatives
			
		mat b = e(b)
		matrix score double xb = b
		
		tab lottype, gen(lt_)
		tab borough, gen(b_)
		tab yearpermit, gen(y_)

		gen pr = exp(xb)/(1+exp(xb))
			
		gen assessprop = pr * ((assesstot - assessland)/(frac_assess*underassess)) * dtaxrate_onsite

		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		
		collapse (sum) assessprop inclusionaryunits
		
		gen acostperunit = assessprop/inclusionaryunits
		
		summ acostperunit
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Marginal fiscal cost

	capture program drop boot_mcost
	program define boot_mcost, rclass
				 		  
		preserve
								
		logit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict xb, xb
		
		gen pr = exp(xb)/(1+exp(xb))
		
		gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
		
		gen mvalue_bldg = (assesstot - assessland)/(frac_assess*underassess)
			
		gen assessprop = pr * mvalue_bldg * dtaxrate_onsite
		gen assessprop_plus = pr_plus * mvalue_bldg * (dtaxrate_onsite + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		gen inclusionaryunits_plus = pr_plus * floor(`sh_inclusionary' * unitsres)
		
		collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)
		
		summ mcostperunit
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Marginal fiscal cost, RE logit

	capture program drop boot_mcost_relogit
	program define boot_mcost_relogit, rclass
				 		  
		preserve
								
		melogit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit || nta_boot:, vce(cl nta_boot) asis
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict eta, eta
		
		gen pr = exp(eta)/(1+exp(eta))
		
		gen pr_plus = exp(eta+`coef_dtaxrate'*0.001)/(1+exp(eta+`coef_dtaxrate'*0.001))
		
		gen mvalue_bldg = (assesstot - assessland)/(frac_assess*underassess)
			
		gen assessprop = pr * mvalue_bldg * dtaxrate_onsite
		gen assessprop_plus = pr_plus * mvalue_bldg * (dtaxrate_onsite + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		gen inclusionaryunits_plus = pr_plus * floor(`sh_inclusionary' * unitsres)
		
		collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)
		
		summ mcostperunit
		return scalar tmp = r(mean)
		
		restore
		
	end

*** Marginal fiscal cost, probit

	capture program drop boot_mcost_probit
	program define boot_mcost_probit, rclass
				 		  
		preserve
								
		probit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict xb, xb
		
		gen pr = exp(xb)/(1+exp(xb))
		
		gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
		
		gen mvalue_bldg = (assesstot - assessland)/(frac_assess*underassess)
			
		gen assessprop = pr * mvalue_bldg * dtaxrate_onsite
		gen assessprop_plus = pr_plus * mvalue_bldg * (dtaxrate_onsite + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		gen inclusionaryunits_plus = pr_plus * floor(`sh_inclusionary' * unitsres)
		
		collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)
		
		summ mcostperunit
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Marginal fiscal cost, LPM

	capture program drop boot_mcost_lpm
	program define boot_mcost_lpm, rclass
				 		  
		preserve
								
		reg inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot)
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict xb, xb
		
		gen pr = xb
		
		gen pr_plus = xb+`coef_dtaxrate'*0.001
		
		gen mvalue_bldg = (assesstot - assessland)/(frac_assess*underassess)
			
		gen assessprop = pr * mvalue_bldg * dtaxrate_onsite
		gen assessprop_plus = pr_plus * mvalue_bldg * (dtaxrate_onsite + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		gen inclusionaryunits_plus = pr_plus * floor(`sh_inclusionary' * unitsres)
		
		collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)
		
		summ mcostperunit
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Marginal fiscal cost, GMM

	capture program drop boot_mcost_gmm
	program define boot_mcost_gmm, rclass
				 		  
		preserve
				
		logit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b_init = e(b)
				
		gmm (inclusionary_onsite - invlogit({xb:dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit constant})), ///
			instruments(dtaxrate_onsite_sim $controls_lot_bs $controls_block i.borough i.yearpermit) vce(cluster nta_boot) ///
			conv_maxiter(10) from(b_init) quickderivatives
			
		mat b = e(b)
		matrix score double xb = b
		
		tab lottype, gen(lt_)
		tab borough, gen(b_)
		tab yearpermit, gen(y_)

		gen pr = exp(xb)/(1+exp(xb))
		gen pr_plus = exp(xb+b[1,1]*0.001)/(1+exp(xb+b[1,1]*0.001))
			
		gen assessprop = pr * ((assesstot - assessland)/(frac_assess*underassess)) * dtaxrate_onsite
		gen assessprop_plus = pr_plus * ((assesstot - assessland)/(frac_assess*underassess)) * (dtaxrate_onsite + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		gen inclusionaryunits_plus = pr_plus * floor(`sh_inclusionary' * unitsres)
		
		collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)
		
		summ mcostperunit
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** All, GMM

	capture program drop boot_all_gmm
	program define boot_all_gmm, rclass
				 		  
		preserve
				
		logit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b_init = e(b)
				
		gmm (inclusionary_onsite - invlogit({xb:dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit constant})), ///
			instruments(dtaxrate_onsite_sim $controls_lot_bs $controls_block i.borough i.yearpermit) vce(cluster nta_boot) conv_maxiter(10) ///
			from(b_init) quickderivatives
			
		mat b = e(b)
		matrix score double xb = b
		
		gen pr = exp(xb)/(1+exp(xb))
		gen pr_plus = exp(xb+b[1,1]*0.001)/(1+exp(xb+b[1,1]*0.001))
			
		gen assessprop = pr * ((assesstot - assessland)/(frac_assess*underassess)) * dtaxrate_onsite
		gen assessprop_plus = pr_plus * ((assesstot - assessland)/(frac_assess*underassess)) * (dtaxrate_onsite + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		gen inclusionaryunits_plus = pr_plus * floor(`sh_inclusionary' * unitsres)
		
		collapse (rawsum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus (mean) pr_plus pr [aw=unitsres]
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)
		
		gen acostperunit = assessprop/inclusionaryunits
		
		gen elast = 1000*(pr_plus - pr)
	
		summ acostperunit
		return scalar acostperunit = r(mean)
		
		summ mcostperunit
		return scalar mcostperunit = r(mean)
		
		summ elast
		return scalar elast = r(mean)
		
		restore
		
	end


*** Marginal fiscal cost, REWEIGHTED

	capture program drop boot_mcost_rw
	program define boot_mcost_rw, rclass
				 		  
		preserve
				
		logit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
				
		predict xb, xb
		gen pr = exp(xb)/(1+exp(xb))
		
		gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
			
		gen assessprop = pr * ((assesstot - assessland)/(frac_assess*underassess)) * dtaxrate_onsite
		gen assessprop_plus = pr_plus * ((assesstot - assessland)/(frac_assess*underassess))  * (dtaxrate_onsite+ 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = pr * floor(`sh_inclusionary' * unitsres)
		gen inclusionaryunits_plus = pr_plus * floor(`sh_inclusionary' * unitsres)
		
		collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus [aw=reweight]
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)
		
		summ mcostperunit
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Elasticity

	capture program drop boot_elasticity
	program define boot_elasticity, rclass
				 		  
		preserve
						
		logit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict xb, xb
			
		gen pr = exp(xb)/(1+exp(xb))
		gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
	
		collapse (mean) pr_plus pr [aw=unitsres]
	
		gen elast = 1000*(pr_plus - pr)
		
		summ elast
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Elasticity, RE logit

	capture program drop boot_elasticity_relogit
	program define boot_elasticity_relogit, rclass
				 		  
		preserve
						
		melogit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit || nta_boot:, vce(cl nta_boot) asis
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict eta, eta
			
		gen pr = exp(eta)/(1+exp(eta))
		gen pr_plus = exp(eta+`coef_dtaxrate'*0.001)/(1+exp(eta+`coef_dtaxrate'*0.001))
	
		collapse (mean) pr_plus pr [aw=unitsres]
	
		gen elast = 1000*(pr_plus - pr)
		
		summ elast
		return scalar tmp = r(mean)
		
		restore
		
	end

*** Elasticity, probit

	capture program drop boot_elasticity_probit
	program define boot_elasticity_probit, rclass
				 		  
		preserve
						
		probit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict xb, xb
			
		gen pr = exp(xb)/(1+exp(xb))
		gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
	
		collapse (mean) pr_plus pr [aw=unitsres]
	
		gen elast = 1000*(pr_plus - pr)
		
		summ elast
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Elasticity, LPM

	capture program drop boot_elasticity_lpm
	program define boot_elasticity_lpm, rclass
				 		  
		preserve
						
		reg inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot)
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		predict xb, xb
			
		gen pr = xb
		gen pr_plus = xb+`coef_dtaxrate'*0.001
	
		collapse (mean) pr_plus pr [aw=unitsres]
	
		gen elast = 1000*(pr_plus - pr)
		
		summ elast
		return scalar tmp = r(mean)
		
		restore
		
	end
	
*** Elasticity, GMM

	capture program drop boot_elasticity_gmm
	program define boot_elasticity_gmm, rclass
				 		  
		preserve
				
		logit inclusionary_onsite dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) asis
		mat b_init = e(b)
				
		eststo: gmm (inclusionary_onsite - invlogit({xb:dtaxrate_onsite $controls_lot_bs $controls_block i.borough i.yearpermit constant})), ///
			instruments(dtaxrate_onsite_sim $controls_lot_bs $controls_block i.borough i.yearpermit) vce(cluster nta_boot) ///
			from(b_init) conv_maxiter(10) quickderivatives
			
		mat b = e(b)
		matrix score double xb = b
		local coef_dtaxrate = b[1,1]
		
		gen pr = exp(xb)/(1+exp(xb))
		gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
	
		collapse (mean) pr_plus pr [aw=unitsres]
	
		gen elast = 1000*(pr_plus - pr)
		
		summ elast
		return scalar tmp = r(mean)
		
		restore
		
	end
	

*** Elasticity, multinomial logit (for offsite/onsite)

	capture program drop boot_elasticity_mnl
	program define boot_elasticity_mnl, rclass
				 		  
		preserve
		
		drop if yearpermit < 2009
				
		capture drop pr1 pr2 xb1 xb2 mult_outcome
		capture constraint drop _all
		
		gen mult_outcome = .
		replace mult_out = 0 if max(inclusionary_offsite,inclusionary_onsite) == 0
		replace mult_out = 1 if inclusionary_offsite == 1
		replace mult_out = 2 if inclusionary_onsite == 1

		* Define constraints from economic theory

		constraint define 1 [#2]:dtaxrate_onsite = 0
		constraint define 2 [#3]:dtaxrate_offsite = 0
		constraint define 3 [#3]:dtaxrate_onsite = [#2]:dtaxrate_offsite
			
		mlogit mult_outcome dtaxrate_onsite dtaxrate_offsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) constraints(1 2 3) iter(30) baseoutcome(0)
		
		mat b = e(b)
		local coef_dtaxrate = b[1,44]
		
		predict pr1, pr outcome(1)
		predict pr2, pr outcome(2)
		
		predict xb1, xb outcome(1)
		predict xb2, xb outcome(2)
		
		gen pr1_plus = exp(xb1+`coef_dtaxrate'*0.001)/(1+exp(xb1+`coef_dtaxrate'*0.001)+exp(xb2))
		gen pr2_plus = exp(xb2+`coef_dtaxrate'*0.001)/(1+exp(xb1)+exp(xb2+`coef_dtaxrate'*0.001))
		
		collapse (mean) pr1 pr1_plus pr2 pr2_plus [aw=unitsres]
	
		gen elast1 = 1000*(pr1_plus - pr1)
		gen elast2 = 1000*(pr2_plus - pr2)
			
		summ elast1
		return scalar elast_offsite = r(mean)
		
		summ elast2
		return scalar elast_onsite = r(mean)
		
		restore
		
	end

*** Marginal fiscal cost, multinomial logit (for offsite/onsite)

	capture program drop boot_mcost_mnl
	program define boot_mcost_mnl, rclass
				 		  
		preserve
		
		drop if yearpermit < 2009
				
		capture drop pr1 pr2 xb1 xb2 mult_outcome
		capture constraint drop _all
		
		gen mult_outcome = .
		replace mult_out = 0 if max(inclusionary_offsite,inclusionary_onsite) == 0
		replace mult_out = 1 if inclusionary_offsite == 1
		replace mult_out = 2 if inclusionary_onsite == 1

		* Define constraints from economic theory

		constraint define 1 [#2]:dtaxrate_onsite = 0
		constraint define 2 [#3]:dtaxrate_offsite = 0
		constraint define 3 [#3]:dtaxrate_onsite = [#2]:dtaxrate_offsite
			
		mlogit mult_outcome dtaxrate_onsite dtaxrate_offsite $controls_lot_bs $controls_block i.borough i.yearpermit, cl(nta_boot) constraints(1 2 3) iter(30) baseoutcome(0)
		
		mat b = e(b)
		local coef_dtaxrate = b[1,44]
		
		predict pr1, pr outcome(1)
		predict pr2, pr outcome(2)
		
		predict xb1, xb outcome(1)
		predict xb2, xb outcome(2)
		
		gen pr1_plus = exp(xb1+`coef_dtaxrate'*0.001)/(1+exp(xb1+`coef_dtaxrate'*0.001))
		gen pr2_plus = exp(xb2+`coef_dtaxrate'*0.001)/(1+exp(xb2+`coef_dtaxrate'*0.001))
		
		gen assessprop1 = pr1 * ((assesstot - assessland)/(frac_assess*underassess)) * dtaxrate_offsite
		gen assessprop1_plus = pr1_plus * ((assesstot - assessland)/(frac_assess*underassess)) * (dtaxrate_offsite + 0.001)
			
		gen assessprop2 = pr2 * ((assesstot - assessland)/(frac_assess*underassess)) * dtaxrate_onsite
		gen assessprop2_plus = pr2_plus * ((assesstot - assessland)/(frac_assess*underassess)) * (dtaxrate_onsite + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits1 = `sh_inclusionary' * pr1 * unitsres 
		gen inclusionaryunits1_plus = `sh_inclusionary' * pr1_plus * unitsres
		gen inclusionaryunits2 = `sh_inclusionary' * pr2 * unitsres 
		gen inclusionaryunits2_plus = `sh_inclusionary' * pr2_plus * unitsres
		
		collapse (sum) assessprop1 assessprop1_plus assessprop2 assessprop2_plus inclusionaryunits1 inclusionaryunits1_plus inclusionaryunits2 inclusionaryunits2_plus
		
		gen mfc_offsite = (assessprop1_plus - assessprop1) / (inclusionaryunits1_plus - inclusionaryunits1)
		gen mfc_onsite = (assessprop2_plus - assessprop2) / (inclusionaryunits2_plus - inclusionaryunits2)
			
		summ mfc_offsite
		return scalar mfc_offsite = r(mean)
		
		summ mfc_onsite
		return scalar mfc_onsite = r(mean)
		
		restore
		
	end
